The Quantum de Laval Nozzle: stability and quantum dynamics of sonic horizons in a 
toroidally trapped Bose gas containing a superflow 
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We study an experimentally realizable system containing stable black hole-white hole acoustic 
horizons in toroidally trapped Bose-Einstein condensates - the quantum de Laval nozzle. We nu- 
merically obtain stationary flow configurations and assess their stability using Bogoliubov theory, 
finding both in hydrodynamic and non-hydrodynamic regimes there exist dynamically unstable re- 
gions associated with the creation of positive and negative energy quasiparticle pairs in analogy with 
the gravitational Hawking effect. The dynamical instability takes the form of a two mode squeezing 
interaction between resonant pairs of Bogoliubov modes. We study the evolution of dynamically 
unstable flows using the truncated Wigner method, which confirms the two mode squeezed state 
picture of the analogue Hawking effect for low winding number. 



I. INTRODUCTION 

The idea of analogue gravity [1], which includes the 
possibility of observing the analogue of the Hawking ef- 
fect in fluid systems exhibiting sonic horizons, was first 
put forward by Unruh [2]. In that paper, using an analy- 
sis similar to Hawking's original analysis for cosmological 
black holes [3, 4], Unruh showed an acoustic black hole 
should emit sound waves with a Planckian spectrum at 
the Hawking temperature 



2'hch 



(1) 



where gn is the surface gravity at the black hole hori- 
zon and ch is the speed of sound at the horizon. The 
derivation required the quantization of a scalar field prop- 
agating in a classical fluid, analogous to a classical gravi- 
tational field. In Bose-Einstein condensates (BECs) long 
wavelength excitations propagate hydrodynamically, giv- 
ing a direct analogue model in a system where the Hawk- 
ing temperature is relatively large [5], and potentially 
measurable using recently proposed schemes based on 
Raman spectroscopy [6] , or parametric resonance [7] . 

Even in the ultra-cold regime where BECs occur, an 
extremely low temperature and low losses are both desir- 
able features of any experimental setup to test analogue 
Hawking radiation in BECs, since it is inherently a deli- 
cate and weak signal. One way to realize a sonic horizon 
in a trapped BEC without any outcoupling is by set- 
ting up a persistent supercurrent in a toroidal trap with 
"bumps" in the trapping potential. The bumps act like 
constrictions, creating a de Laval nozzle [8] type configu- 
ration (Fig. 1). The experimental realization of a toroidal 



magnetic trap for ultracold atoms, first demonstrated by 
using two current-carrying loops [9], has since been de- 
veloped using magnetic waveguides [10], a microchip trap 
[11] and a four loop configuration [12]. The requisite 
bumps could be introduced optically by a detuned laser 
shining through an appropriately patterned mask [13]. 
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FIG. 1: Hydrodynamic de Laval nozzle. When the flow 
speed (v) of a normal fluid through a constriction achieves 
the sound velocity (c) at the waist, the outgoing flow becomes 
supersonic. To maintain continuity, fluid density (grayscale) 
and stream line spacing (solid lines) decrease from left to right 
as the fluid accelerates. 

While there exist previous theoretical studies of de 
Laval nozzle geometries in the context of acoustic black 
holes, both for classical fluids [14, 15] and for BECs 
[5, 16], the analyses in these studies are either classi- 
cal or semiclassical in nature. Previous investigations of 
acoustic black hole geometries [5, 16, 17] have focused 
on regimes where both the hydrodynamical and geomet- 
ric acoustics descriptions for a BEC apply, so that the 
semiclassical WKB method can be applied to calculate 
the Hawking temperature in close analogy with the grav- 
itational derivation [2-4, 18]. In particular, flows are 
treated as hydrodynamic, and the effects of the trap are 
neglected through some form of local density or WKB 
approximation. This is understandable given the con- 
ditions under which Hawking first discovered the effect, 
but our primary interest are quantum effects which have 
also been studied for other BEC acoustic horizon scenar- 
ios [5, 16, 17, 19-21]. 

In this work, we consider a quite different regime where 
hydrodynamics is valid only for the low energy modes, 
whereas geometric acoustics is only valid in the limit of 
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high winding number. The system of interest is in the 
region of parameter space where hydrodynamics and geo- 
metric acoustics approximations are not usually applica- 
ble, but where progress can be made with more detailed 
numerical analysis. 

We introduce and analyze a system which exhibits an 
acoustic black hole and white hole horizon in a trapped 
BEC, formed by two de Laval nozzles in a toroidal ge- 
ometry. Under the conditions of steady flow this config- 
uration represents the simplest stable de Laval geometry 
for a Hamiltonian BEC system, and has some appealing 
properties for studying the analogue Hawking effect. In 
particular, it has a discrete excitation spectrum and peri- 
odic boundary conditions. Our primary aims are to find 
stationary solutions of the Gross-Pitaevskii equation for 
this system, and to investigate their stability and quan- 
tum dynamics. 



II. QUANTUM DE LAVAL NOZZLE 

A weakly interacting Bose gas trapped in a one di- 
mensional potential V{x) at zero temperature is well de- 
scribed by the Gross-Pitaevskii equation [19, 22] 
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where the effective one dimensional interaction strength 
is Ui-D = UQ/{'iTTr\) and Uq = 'i-Kh^as/m is the usual 
s-wave interaction parameter. The reduction to one di- 
mension assumes the transverse wavefunction is in the 
harmonic oscillator ground state of toroidal trap: (j){r) = 



,1/2 
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where ujj_ is the transverse trapping frequency. The dy- 
namics become effectively one dimensional in this regime 
since the transverse motion is frozen out by the large 
energy required to excite transverse modes. The s-wave 
scattering description remains valid provided the scatter- 
ing length a is much smaller than the transverse dimen- 
sion r_L so that the scattering remains effectively three 
dimensional [23]. The wavefunction can be written as a 
macroscopic order parameter 



with current density 
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and velocity v = hdxQ/rn. In density-phase variables 
the system is governed by the equations of motion 
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and 
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When the interaction term dominates, the density varies 
slowly, and the Laplacian term (quantum pressure) 
can be dropped; then the last equation can be takes the 
form of Euler's equation 
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and we recover the classical isentropic flow equations. 
Combining (5) and (7), one can derive the nozzle equation 
[8, 24]. 
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relating variations of the potential and the flow veloc- 
ity. The physical consequences of this form of the noz- 
zle equation are as follows (refer to Fig. 1): Sonic flow 
(v = c) is only permitted where dV = 0, that is at 
the waist of the nozzle. For subsonic flow (v < c), 
when dV is negative/positive the velocity is decreas- 
ing/increasing. Conversely for supersonic flow, when dV 
is negative/positive the velocity is increasing/decreasing. 
Therefore, if the flow is subsonic on approach to the noz- 
zle waist, becoming sonic at the waist, it becomes su- 
personic on exiting the waist, and conversely for an ap- 
proaching supersonic flow. 

To achieve transonic steady flow in a toroidal geometry, 
two de Laval nozzles are required in tandem, the flow 
becoming supersonic at the waist of the first, and then 
subsonic at the waist of the second. This configuration 
corresponds to the formation of both a black and white 
hole horizon. To implement such a geometry we consider 
an external potential of the form 
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which has periodicity 2 over the interval —L/2 < x < 
L/2, which is periodic for the toroidal geometry. For 
a BEC confined by such a potential, the stationary 
states are found by solving the time-independent Gross- 
Pitaevskii equation subject to phase quantization. Since, 
as we will see below, the solutions exhibit strong modifi- 
cations to hydrodynamic behavior, we refer to this con- 
figuration as the quantum de Laval nozzle (QdLN). 

We wish to study currently realistic or potentially 
achievable parameters. Toroidal ultra-cold atom waveg- 
uides have been developed by several groups [10, 12]. 
We take as nominal values those of the recent experi- 
ments of the Stamper-Kurn group [10]. The experiments 
typically consist of iVo '--^ 3 x 10'"' ^^Rb atoms held in a 
toroidal trap with transverse frequency luj_ ~ 2tt 80 Hz, 
and radius i? ~ 1 mm which gives an azimuthal length 
L ^ 6.2 mm and a transverse harmonic oscillator dimen- 
sion r± = y^huj±/m ^ 1.2 fj.ni. Typical winding num- 
bers are estimated by assuming that the circulation veloc- 
ity is constant and assuming that the gas fills the entire 
perimeter of the toroid (which is not the case in the ex- 
periment). This leads to the estimate wq = mLp' /2iThT, 
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which for the circulation period of Ref. [10] 200ms) 
gives Wo ~ 3 X lO** which is quite large. Construct- 
ing stationary solutions for the quantum de Laval noz- 
zle for these exact parameters presents a major com- 
putational challenge as it is necessary to resolve phase 
variations of the wavefuction on a very small scale, how- 
ever we note that the condensates in this experiment are 
launched into the toroid and do not form periodic sta- 
tionary solutions. In theoretical work existing in the lit- 
erature more modest winding numbers have been consid- 
ered: 1 < Wo ^ 10 [19]. In this work we will examine the 
stability of a similarly modest range of winding numbers 
in detail, and also for comparison we include wq = 50. 



III. STATIONARY STATES 

In this section we find stationary solutions for the 
QdLN using the Gross-Pitaevskii equation and com- 
pare the results with hydrodynamic and perturbative ap- 
proaches. 

It is convenient to normalize the condensate wavefunc- 
tion to the single particle form J |-0(a;)|^ da; = 1 in what 
follows, so that g = NoUir, hereafter describes the to- 
tal effective nonlinearity for A^o condensate atoms. For 
steady state fiow wc take the one dimensional stationary 
solution of the form 



(10) 



If wc take the stationary solution of (5) and (6) using 
(10), and take the fixed current condition, wc can write 
n = s'^,v = J/s'^, 



(11) 



Solutions to this nonlinear equation then allow us to re- 
construct the wavcfunction by specifying: 

J dy 



'ip{x) — s(x) exp («7?(x)) . 
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Solutions must also satisfy the phase quantization condi- 
tion 

f.L/2 

v{x)dx = 27ru;o (14) 
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for an integer winding number given by wo- 
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A. Hydrodynamic solutions 

When the interactions dominate the density varies 
slowly so we can invoke the hydrodynamic approxima- 
tion and drop the Laplacian term. In this case (11) can 
be written as a cubic, either in terms of the density: 
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or in terms of the velocity: 
+ 2' ^ ^ ' 
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For the case where the flow is zero, there is one non- 
trivial solution to (15) 



(17) 



The chemical potential is /i = ^/L + Vo/2, which yields 
the density 
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On the other hand, for the case where there is non- 
zero flow (J > 0), we find solutions using (16) since the 
equations have a simpler form in this case. The solu- 
tions are conveniently separated by the discriminant of 
the cubic [25] 
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Transonic configurations exist when there are two real 
positive solutions, which occurs when d(x) < for all x. 
Since {Jg/m)^ > the negative semi-definite character 
of d{x) imposes the constraint Vq < fj,. We can express 
these solutions analytically as: 



with 
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(20) 
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Note v^{x) is the subsonic branch, whereas v+{x) is the 
supersonic branch. 

A continuous, single valued transonic solution can be 
constructed when the two positive solutions coincide at 
the horizon, which we take to be a; = xh = 0. This oc- 
curs when d[xH) = 0. At the horizon V{xh) = (ie. 
the maximum of the potential acts as the waist of the 
de Laval nozzle), so rearranging (19) we find the condi- 
tion for transonic flow is given by the critical chemical 
potential 



A^crit 



,1/3 



(23) 



Note for /i < ^crit we have d{x) > and the flow is 
unstable. 

We can find a transonic solution by taking the chem- 
ical potential ^ = ficntixH) so that there is a crossover 
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from the subsonic to supersonic branches at the horizon: 
V-{xh) = Vj^(xh)- The transonic solution is then con- 
structed by joining the subsonic [v-) and supersonic {v^) 
solution branches. Without loss of generality, we take the 
subsonic branch to span the interval x € [—L/2, 0], with 
the supersonic branch in the interval x £ [0,L/2]. Consis- 
tent solutions with integer winding number wq are found 
by iterating the hybrid solutions and the constraints to 
find the appropriate conserved current J n{x)v{x). 
The resulting stationary solution is fully determined by 
the parameters Vq, g and Wq. 



B. Perturbation theory 

Although we have the analytical solutions of the hy- 
drodynamic theory it is useful to adopt a perturbative 
approach which has the advantage of giving simple and 
reasonably accurate solutions at first order in powers of 
fi/2 ^ (Vb/Ato)^^^, where /ip is the zeroth order chemi- 
cal potential. We have chosen our potential so that we 
can choose the unperturbed problem as the homogeneous 
solution for Vb = 0, with critical flow so that v = c ev- 
erywhere. 

At any order the solutions must satisfy the set of equa- 
tions 

2 ^ . n f2'Kx\\ 2Jg 



0, (24) 
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At zeroth order, the potential free equations satisfied by 
the unperturbed variables (/ig, Jojfo) are 

3 2fiQVo , 2Jo5 
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In solving the cubic we find solutions vq € 
(-2( Jog/m)i/3, {Jog/my^^ ( Jo5/m)i/3). Only the pos- 
itive flow solutions are physical and their coalescence at 
zeroth order is helpful at higher order where the solutions 
break the parity symmetry of the potential. 

At zeroth order the solutions can expressed in terms 
of the winding number as fo — 2nhwo/mL, Jq = rnvg/g, 
/io = 3touq/2, and uq ~ niv'^/g. Since we are going to 
require Vq <C /io this imposes the condition 

'^y^L^\"\, ton, 



for the validity of the perturbation series. We introduce 
the rescaling (w, J, a;) = (vvq, JJq, flfj.o,xL), to obtain 
the equations 
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The repeated solution at zeroth order means we have to 
use a perturbation series in powers of e^/^ [26], so we as- 
sume an expansion of the form v = 1 + e^^^vi 



ev2. ■ ■ 

and similarly for the other variables. We can obtain con- 
sistent solutions to Eqs. (33-36) up to 0(e) which give 
a good qualitative description of the solutions and are 
quite accurate for a wide range of parameters. 

Terms in the expansion of (33) of order e°, e^/^ cancel, 
and the e equation is 



/ii + sin^ (27rx). 



(37) 



Substituting the series into (34) gives yUi = which is 
not surprising, and in fact fj,2 = 0- The subsonic and 
supersonic solutions are automatically matched at the 
acoustic horizon by choosing vi ~ sin {2ttx), whereby the 
supersonic region is < a; < 1/2, coinciding with our 
previous hydrodynamic treatment. 

Returning to dimensioned variables, the first order so- 
lutions for the quantum de Laval nozzle are 



. = .o|l+|^j sin(^ 



n = uq 1 



Vo \ . / 2Trx 
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with J = Jo+OiiVo/fiof/^) and ^l = Mo+O((yo/Mo)'/')- 
These expressions give a good qualitative description of 
the Gross-Pitaevskii solutions, and are typically very 
close to the full hydrodynamic solutions (see Fig. 2). As 
expected, the differences are more apparent with increas- 
ing e corresponding to a more important quantum pres- 
sure term, and with increasing distance from the sonic 
horizons. 



C. Solutions of the Gross-Pitaevskii equation 

We now use the transonic solution to the hydrody- 
namic problem (Sec. Ill A) as a starting point to find- 
ing the stationary solutions of the Gross-Pitaevskii equa- 
tion (GPE). The full numerical solutions exhibit "ripple" 
structures in regions where the quantum pressure term 
becomes important, in particular, in the region of super- 
sonic flow downstream of the black hole acoustic horizon. 

In order to find stationary solutions we use constrained 
optimization. Imaginary time evolution is not feasible in 



5 



1.5 - 



CM 



0.5 





GPE 

perturbation 

hydrodynamic 


'/ 


(a) 



3 




FIG. 2: (color online) QdLN stationary flow. Parameters 
are wo = 3, Vo = lOOftcji (energies are in units of hujL = 
h^/mL^), fi = 5.93 X lO^fttJi and g = 3.51 x lO^LRtJL. (a) 
Comparison of solutions for condensate density, (b) Velocity 
and speed of sound for GPE solution. 



this case because we are interested in stationary states 
which are excited into circular motion relative to the 
ground state. We formulate the problem as the the min- 
imization of the Gross-Pitaevski functional for a fixed 
nonlinearity g, potential depth Vq and current J, subject 
to a phase quantization constraint in terms of a fixed 
winding number wq- The problem is recast as the set of 
algebraic equations: 



7^-r^+ tJ'-y{x)- gs[x) - ^ 



m 



J 



S(x)2 



six) = 0, 
dx — 27rw = 0. 



The phase circulation constraint ensures the wavefunc- 
tion ■(/; = s(x)e'''^^-' is everywhere single- valued. 

The solution for our vector of unknowns X = {s{xi), /i} 
is found by Lcvenberg-Marquardt optimization [27, 28] in 
MATLAB using the hydrodynamic solution as the initial 
condition Xg. The unit of energy for this system which 
we will use to display our results is Hlvl = li' /{mL'^). 

To consider some examples we use a potential with 
Vb = 100?iaJL, and show solutions for winding numbers 
Wo = 3 (Fig. 2) and ~ 10 (Fig. 3). In each case, the 
full hydrodynamic, first order perturbation theory and 
GPE solutions are shown, as well as the flow velocity 
and speed of sound for the GPE stationary solution, and 
the location of the acoustic black hole (BH) and white 
hole (WH) horizons. The ergoregion [v > c) is given 
approximately by the right hand region < a; < 0.5. 
By construction, in the hydrodynamic case, the black 
hole horizon occurs at xbh = 0, whereas the white hole 
horizon occurs at .twh = ±0.5. 

In principle we might expect to be able to vary the 
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FIG. 3: (color online) QdLN stationary flow. Parameters 
are as in Fig. 2, but with = 10, ^ = 5.99 x IO^Twjl and 
g = 3.95 X lO^LftwL. 



winding number wo and potential depth Vo to find a con- 
tinuous range of transonic solutions for the QdLN. This 
is the case, for example, for the toroidal system consid- 
ered by Garay et al. [19], where for a given wq a stability 
diagram over a continuous range of Vq and g was mapped 
out. In fact we find that for the QdLN the total nonlin- 
earity g is not a free parameter of the solutions - it is 
uniquely determined for a given (wo,Vb). This is phys- 
ically reasonable because wq sets the flow velocity and 
the nonlinearity determines the speed of sound, and the 
two must be equal at the sonic horizons. 

We can easily flnd a very accurate relation between 
nonlinearity and winding number using the zeroth or- 
der perturbation theory: g ~ LhujL{2TTWo)'^ . We should 
expect deviations from this relationship for low winding 
number due to the importance of the quantum pressure 
term, however we flnd they are very small. Numerically 
wc flnd for the GPE solutions that g varies according to 
this quadratic law and depends only very weakly on Vq. 
For the values of (wq, Vq) used in this paper the behavior 
is essentially independent of Vq and shows a maximum 
deviation from the zeroth order perturbation theory re- 
sult of order 1%, occurring at our lowest winding number, 
Wq = 3. 

Qualitatively, the main departure of the hydrodynamic 
solutions from flrst order perturbation theory is a loss 
of parity: the increasingly important interaction energy 
eventually lifts the antisymmetry of vi ~ sin (27ra:;/L); 
such differences are more pronounced for low winding 
number. 



IV. QUASIPARTICLES AND STABILITY 

We will now determine the stability of the GPE sta- 
tionary solutions by flnding their Bogoliubov excitation 
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spectra for a wide range of potential depths and winding 
numbers. For unstable configurations the standard Bo- 
goliubov analysis is insufficient, and we use the theory of 
Leonhardt et al. [17] to obtain the correct normalizable 
Bogoliubov modes. These modes show some interesting 
localization properties with respect to the acoustic hori- 
zons, and will be used to construct Bogoliubov vacuum 
states when we come to dynamical simulations in Sec. V. 
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A. Normalizable Bogoliubov modes 
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The linear excitations of the condensate are described 
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by the Bogoliubov-de Gennes (BdG) equations [29, 30]. 
Consider a solution with small oscillations around a sta- 
tionary state 
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where f3i and /3* are the amplitudes for the oscillations 
and 4'oix) is normalized to unity in what follows. For 
the quantum field, these quantities are replaced by the 
bosonic annihilation and creation operators for the exci- 
tations, given by bi and b] respectively, with [hi, Sj] — Sij. 
4>o{x) is the solution to the time- independent GPE given 
by (2). To obtain the correct physical modes it is neces- 
sary to introduce the operator 



Q=l-|0o)(0o| 



(41) 



which projects orthogonal to the condensate. Introduc- 
ing projectors appropriately and substituting the Bogoli- 
ubov expansion into the time-dependent GPE and while 
only keeping terms linear in Ui{x) and Vi{x) yields the 
modified BdG equations 



(42) 
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where the operator £ is given by 
£GP-Ai+.9QI0oPQ 



£ = 



£gp-M- 



and the Gross-Pitacvskii operator is 



£gp 



2m 



Vix)+g\M 



(43) 



The stationary solution satisfies (£gp — m)0o = 0. The 
solutions to this equation are the eigenvalues = h{uji + 
i"fi), and the normal modes of the system. 

The orthogonality and symmetry relations are fixed by 
the requirement that the many body Hamiltonian for the 
interacting Bose gas is diagonal (to quadratic order) in 
quasiparticle operators and that the transformation to 
quasiparticles preserves the commutation relations. 



The operator £ is not Hcrmitian so that complex eigen- 
values are allowed, corresponding to dynamical instabili- 
ties of the system. It is straightforward to show that the 
modes with complex eigenvalues have zero norm [31] and 
cannot therefore be associated with bosonic operators in 
the field expansion (40). However, following Leonardt 
et al. [17], it is still possible to construct normalizable 
modes for the unstable modes by the construction: 



(44) 
(45) 



where [uf {x) , (x)] and [u~ (x) , v~ (x)] are the eigenvec- 
tors associated with the unstable positive and negative 
energy eigenvalues respectively. Note that due to the 
symmetries of (43) we have uj~ = —ijj'l- Hereafter we 
use Ui{x) and Vi{x) for the full set of orthonormal modes. 
The quadratic Hamiltonian for the stable modes takes the 
standard form for independent harmonic oscillators. The 
creation and annihilation operators for the new modes 
satisfy the commutation relations for bosonic operators, 
but show up as non-diagonal terms in the Hamiltonian 
subspace for the dynamically unstable modes as [17] 

H2 = Y.hu,Mb]j>,+-b]_b,^) 



^ i h-ij \{bj+bj--b]J)]^) 

dxiu;v--u;*vr) 



(46) 



where the sum is taken over only the dynamically un- 
stable modes, and where bj± is the annihilation operator 
and Sjj. the creation operator corresponding to the nor- 
malizable modes (44) and (45). For stable modes the 
hamiltonian reduces to the usual diagonal Bogoliubov 
form H2 = ^^jiblh + 1/2- J dx \Vj\^). Dynami- 
cally unstable modes are therefore associated with non- 
degenerate parametric amplification [32], which leads to 
growth in the unstable modes at the expense of the con- 
densate mode. For short time dynamics, the complex 
eigenvalue will generate exponential growth in each un- 
stable mode. It is this effect that has been suggested 
to provide the closest analogy with the Hawking effect 
for BECs [17, 19]. However, this picture neglects higher 
order interactions that may be present in the full Hamil- 
tonian, and therefore is likely to fail for dynamics on long 
time scales. We will investigate this further in Sec. V. 



B. Stability and mode structure 

For a dynamically unstable configuration, we con- 
struct normalizable modes using the procedure outlined 
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FIG. 4: Eigenspectrum: We show — h{uji + i'yi) for wo — 
3. Circles are for the left axes, and modes are numbered in 
ordered of increasing Ui. Mode 4 is the first positive frequency 
mode, (a) Stable state for Vb = 140?icJi,. (b) Dynamically 
unstable state for Vb — 163.7?ttJi,; the dashed line denotes the 
mode pair with uji — —uje signaling the dynamical instability. 

in Sec. IV A. We additionally sort the eigenvalues in as- 
cending order by and label the modes accordingly. 

Figure 4 shows the eigenvalue spectrum for the first few 
modes with wq — 3 and two different values of Vq: (a) 
Vo = UOhuj^; and (b) Vb = 163.7?iw^. While both cases 
have negative eigenvalues, indicating energetic (Landau) 
instabilities due to non-zero flow, only case (b) exhibits 
dynamical instabilities also. Following the theory of the 
previous section, the onset of a dynamical instability is 
associated with a pair of modes (labelled by j and k 
say) with complex eigenvalues that satisfy ujj ~ — i^fc. 
In particular, case (b) indicates that modes 1 and 6 are 
unstable with uii = —luq. 

Figure 5 shows the stability diagram for the QdLN that 
results by performing the diagonalization for a range of 
parameters, Vq and wq. For each point we have calculated 
the maximum of the absolute value for the imaginary part 
of all eigenvalues. The essential features we observe are: 
(i) there are regions exhibiting dynamic instabilities; (ii) 
these regions become narrower and smaller in magnitude, 
but more closely spaced for larger values of the winding 
number wq, whereas they become broader and larger in 
magnitude as the potential depth Vq increases. 

In Fig. 6 we show the corresponding mode func- 
tions Ui{x) and Vi{x) that result from the solutions of 
the BdG equations for the parameters wq = 10 and 
Vq ~ 139.2fkoL, which has a dynamical instability for 
the modes i = 5, 6. We note, although not shown here, 
the mode functions for wq = 10 and Vb = lOOfiujj^ (dy- 
namically stable) are very similar to the Vb = 139.2tiLUj^ 
case. 

The mode functions exhibit a rich structure for the 
low energy modes, not least being the sort of "local- 
ization" of modes that is associated with the acoustic 
black hole geometries. In particular, for modes i < 4, 
Ui{x) is localized in the region < x < 0.5, which corre- 
sponds to the supersonic region, whereas Vi{x) is local- 
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FIG. 5: Stability diagram: We plot max{|7, |} against Vb 
for stationary solutions with toroidal flow over a range of 
winding numbers. The quantum de Laval nozzle is unsta- 
ble at the narrow spikes in the imaginary eigenvalues (which 
are equally narrow on a logarithmic scale); elsewhere the flow 
is stable. 



ized in the region —0.5 < x < 0, corresponding to the 
subsonic region. For modes i > 7 we find the reverse 
is true in general, although the localization occurs to a 
lesser extent. Modes z = 5, 6 (the dynamically unstable 
modes) indicate a crossover between these two regimes. 
In Fig. 7 we show the average position of the quasiparticle 
modes for two unstable cases for comparison. The mean 
quasiparticle position for each mode is calculated as [33] 
(x)n = {{x)u„ + {x)v„ - {x)^o)/ 1 dx |u„P -l- |u„p, wherc 
{x)un ~ j dx u!^{x)xun{x) with (x)„„ defined similarly, 
and (x)^o = J dx (l)o{x)* X(j)o{x) / J dx \(j)o{x)\'^ . We note 
that the negative energy modes (1-3 for wq = 3 and 1-5 
for Wq ~ 10) are always located significantly inside the 
sonic horizon (x) > 0, while higher energy modes become 
located nearer the horizon. 
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FIG. 6: (Color online) Orthogonal mode functions for dy- 
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the solid blue curve, and lm{Ui{x)) is given by the dashed 
red curve. In column (c) 'Re{V*{x)) is given by the solid blue 
curve, and lm{V*{x)) is given by the dashed red curve. 
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V. DYNAMICS 

To investigate the time dynamics of the QdLN we will 
use the stationary states as our starting point for quan- 
tum field theory simulations using the truncated Wigner 
method. We will compare the quasiparticlc population 
dynamics in both stable and unstable regimes. First we 
briefly discuss the connection between our Bogoliubov 
analysis and the analogue Hawking effect. 



A. Two-mode Bogoliubov model 

The interaction Hamiltonian for a pair of dynamically 
unstable modes with uji ~ —^j, and with |7i| = |7j| = 7 
is 



TJint = ifil[b-b+ - blbl] 



(47) 



which describes the formation of a two-mode squeezed 
state. By finding the time evolution for the two mode 
density matrix according to i/jnt and averaging over the 
negative energy mode we obtain the density operator for 
the positive energy mode 



P+{t) 



1 



cosh 



^(tanh2 7t)" \n) 



(48) 



a thermal state with mean occupation = sinh^ 7i. 
We thus have a loose analogy with the Hawking effect: 
pairs of quasiparticles can be produced with no energy 
cost, such that one quasiparticle enters the negative en- 
ergy state which is located inside the supersonic region 
(for our = 10 case in Fig. (6) this is mode 1). and 
the other is promoted to positive energy (mode 6), which 
is centered much closer to the horizon. Tracing over the 
negative partner gives a thermal state for the postive en- 
ergy mode. This connection has been pointed out previ- 
ously [17], but to our knowledge it has not been confirmed 
for the trapped Bosc gas using analysis of the GPE as we 
are able to do here (see Fig. 9). However, we note that 
the analogy with the gravitational Hawking effect is in- 
complete, as may be seen from the time dynamics of the 
positive energy mode occupation number. A more direct 
Hawking analogue would generate a time independent 
multimodc thermal emission spectrum, whereas here we 
have an exponentially growing, single mode emission. 

Thus far, wc have found the elementary excitations for 
the QdLN and found that this indicates dynamically un- 
stable configurations for certain sets of parameters. In 
order to verify that such configurations do indeed lead to 
exponential growth in the unstable modes, we consider 
the dynamics of the system. To do this we use the trun- 
cated Wigner method to perform short time simulations 
of the full interacting quantum field theory describing the 
trapped Bose gas. 
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B. Truncated Wigner method 

The Wigner representation provides a symmetrically 
ordered formalism for phase space simulations of quan- 
tum field theory. Symmetrically ordered operator aver- 
ages are computed by ensemble averaging many classical 
field trajectories. The truncated Wigner method [34-40] 
involves neglecting intractable third order derivatives in 
the equation of motion for the Wigner distribution. The 
method then reduces to numerically evolving a multi- 
mode classical field using the GPE (2) [34]. The theory 
differs from pure mean field theory in that statistical fluc- 
tuations in the initial state reproduce quantum fluctua- 
tions in the observables extracted by ensemble averaging. 
The method is known to be accurate for short evolution 
times [34]. In the low temperature regime ksT e^, the 
initial field is given by 

,Pix, t = 0) = M^) + E iU^ix)/3, + V*ix)p:] (49) 

where ipQ (x) is a stationary state of the GPE (for our pur- 
poses, the transonic solutions of the QdLN), and where 
Ui and Vi are the Bogoliubov mode amplitudes of the sys- 
tem. The complex random variables Pi arc constructed 
as Pi{t = 0) = {t]i + ir]2)/'/2, where 77, are real, normal 
Gaussian variates with 777 = and rjiffj = Sijifii + 1/2), 
and fii = (e'^i/^sT _ j^^-i jg j-j^g thermal quasiparticle 
occupation. The notation rj represents the stochastic av- 
erage over many samples of ry. In this work we restrict 
our attention to the zero temperature case to investigate 
the stability and dyanamics of our stationary solutions 
in the presence of vacuum fluctuations. 

We expect the details of the quantum dynamics to de- 
pend sensitively on any instabilities, and indeed, accord- 
ing to the two mode model, instabilities can generate 
squeezing. We use the quasiparticle occupation num- 
bers to look for confirmation of this effect in our sim- 
ulations. Dynamically, the Bogoliubov amplitudes can 
be extracted from the classical field as 

P,{t) = J dx{U:{x)i;{x,t) - V:{x)r{x,t)) (so) 

which we use to monitor the populations during our sim- 
ulations. The quasiparticle number in each mode is then 

N,{t)^{brh) = p*mit)-^ (51) 

where the bar indicates an ensemble average over many 
Wigner trajectories. 

In any simulations we must use a restricted basis, and 
our GPE evolution is numerically projected at each time 
step to ensure that the system remains in the low en- 
ergy subspace determined by our energy cut-off. For- 
mally we are using the projected GPE (PGPE) [41] 
to ensure consistent evolution of our restricted phase 
space. For the evolution of the PGPE we have used 



the fourth-order Rungc-Kutta in the interaction pic- 
ture (RK4IP) algorithm [42], adapted to project into 
the low energy subspace defined by our momentum cut- 
off at Ec = h^K^/2m [43]. For the simulations pre- 
sented here we have use Nq = 10^ condensate atoms, 
and M = 1024 modes for the system in the low energy 
subspace, corresponding to a dimcnsionlcss momentum 
cutoff Kc = 2'kM/L = 6.4 x 10'^. For all simulations wc 
have used a time step for the RK4IP algorithm ensur- 
ing the change in total field normalization during each 
trajectory was AN/N < 10"^. 

C. Results 




mode 



FIG. 8: Bogoliubov mode populations for single trajectory 
from TWA evolution for parameters Vb = 163.7hujL, wq — 3, 
corresponding to a dynamically unstable configuration. 
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FIG. 9; Population of the positive energy mode (mode 6 in 
Fig. 10, solid line) with the result from Bogoliubov theory 
{n)+ = sinh^ 'j^t (dashed line). 

For a winding number of = 3, we have carried out 
time dynamical simulations for 40 trajectories using the 
truncated Wigner method. The ensemble averaged quasi- 
particle mode occupations have been calculated for two 
cases (refer to the stability diagram in Fig. 5): (i) For 
the stable case Vq = lAOhcUj^, the quasiparticle modes re- 
main unnocupicd during the interval < w^t < 1. This 
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FIG. 10: Bogoliubov mode populations for single trajectory 
from TWA evolution for parameters Vo — WS.JhcjL, wo — 3, 
corresponding to a dynamically unstable configuration. 



confirms the stability of the system to the extent possible 
given that our GPE stationary state with vacuum noise 
is only an approximation to the true many body sta- 
tionary state, (ii) Fig. 8 shows the populations for the 
unstable case, Vq = IdS.ThuJi,. In the latter case, modes 
1 = 1,6 are unstable, and we observe exponential growth 
in these modes. The growth is seeded by the quantum 
vacuum fluctuations in the initial state and confirms the 
expectation that the system obeys the Hamiltonian (46) 
for a non-degenerate parametric amplifier at short times. 
Note that the mean quasiparticle occupation for modes 
I > 10 is also negligible compared with i < 10. In Fig. 9 
we plot the population in the positive energy mode of the 
unstable pair, iVg, which is seen to be in close agreement 
with the dynamics expected from the Bogoliubov theory 
of Sec. IV. 

We have also investigated the behaviour of the dynam- 
ically unstable configuration for longer times. Single tra- 
jectory results for a simulation time of w^t = 5 are shown 
in Fig. 10. Here we observe growth in the unstable modes 
(i = 1, 6) until lulI 3.25 where there is a peak in 
the mode populations, followed by a decay of occupation 
numbers. Therefore the system undergoes a period of ex- 
citation followed by an apparent return to the initial un- 
excitcd state. This is also evident in the coordinate space 
density plots for the same simulation given in Figs. 11 (a) 
and (b). In particular, plot (a) shows large scale density 
fluctuations for 2 < t < 3.5. Plot (b) shows the den- 
sity relative to the initial state, from which it is clear 
that the density fluctuations are localised in the region 
< a; < 0.5 corresponding to the supersonic region for 
the system. It seems plausible that such excitations and 
revivals should continue to repeat, which would result in 
a "ringing" type excitation of the condensate. However, 
due to the significant computational time required, we 
did not check this prediction. 

The recurrence of the system is evidently due to nonlin- 




FIG. 11: (Color online) Coordinate space density vs time 
for unstable configuration with wq = 3 and Vo = 163.7hu)L'- 
(a) gives the normalized density of the field by n{x, t) = 
\jp{x,t)\'^ /No; (b) gives an intensity plot of the change in den- 
sity from the initial state by An{x,t) = n{x,t) — n{x,Q). In 
plot (a) the high frequency noise has been filtered out for 
clarity. 



ear mode mixing, which is neglected in the BdG analysis. 
In particular, the back-reaction of quasiparticle modes on 
the condensate should become significant for large mode 
occupations. Moreover, the topological constraint im- 
posed by the periodicity of the system (ie. by a fixed 
winding number) means that the decay of circulation 
for the superfluid flow is forbidden. Evidently the de- 
cay channel for this instability is inhibited. The system 
cannot reach a quasi-stationary state corresponding to a 
different value of wq, without a corresponding change in 
Vb and the damping of topological charge; one mecha- 
nism for this process would be soliton shedding, but we 
have not observed this in our simulations of this system. 

We also investigated the dynamics for the higher wind- 
ing number wq ~ 10. In particular, we have examined 
two flows, which are (i) the dynamically stable (accord- 
ing to the BdG analysis) flow, Vq = lOQTiujj^; and (ii) the 
dynamically unstable flow Vq ~ 139.2?iu!j^. We briefly 
summarize the results for this case as we have found it 
to be typical of high winding number behaviour. From 
our ensemble averaged time dynamics carried out in a 
similar manner to the previous cases, we found (1) that 
linear stability, evidenced by time-independent vacuum 
quasiparticle population, was confirmed for short time 
quantum dynamics; (2) for the unstable configuration, 
chaotic multimode dynamics is evident in the otherwise 
stable time interval, in particular, we did not observe 
simple two mode squeezing dynamics of the kind seen 
for wq = 3. This behaviour is not unexpected since for 
the higher winding number, the effective nonlinearity re- 
quired is very large (recall from the perturbation theory 
of Sec. Ill the nonlinearity scales with the square of wind- 
ing number). In our quantum dynamical simulations ele- 
mentary excitations interact with each other, giving rise 
to Landau-Bcliaev damping [44, 45], and this effect is 
more pronounced at higher nonlinearities. 



11 



VI. CONCLUSIONS 

We have introduced and analyzed the quantum do 
Laval nozzle, a toroidal geometry for a BEC that ex- 
hibits both a black and white sonic horizon. Using hy- 
drodynamic theory we have found transonic solutions, 
which we used to find transonic stationary solutions of 
the Gross-Pitaevskii equation. The qualitative proper- 
ties of the GPE solutions are well described by hydrody- 
namic perturbation theory at lowest order. The system 
has broad dynamical instabilities for certain values of the 
winding number wo and potential depth Vq. 

We constructed normalizable Bogoliubov modes for 
the dynamical instabilities, which couple modes of pos- 
itive and negative energy. This analysis leads to a two 
mode squeezing Hamiltonian term corresponding to non- 
degenerate parametric amplification, which leads to ex- 
ponential population growth of unstable mode occupa- 
tion with time - this represents the closest analogy with 
the Hawking effect for our trapped quantum system. 

To analyze this picture further, we have investigated 
the dynamics of several configurations using the trun- 
cated Wigner method which from an analogue model 
point of view includes the effects of nonlinear interac- 
tions between modes and back reaction. For low winding 
number we observe non-degenerate parametric amplifi- 
cation type dynamics at the instability, confirming the 
two-mode Bogoliubov model validity, while for the sta- 



ble configuration there is negligible growth in all modes. 

From the stability analysis, wc note for large wind- 
ing number solutions: (i) the number of unstable regions 
increases, and they become narrower; (ii) the nonlincar- 
ity increases so that the system approaches the hydrody- 
namic regime; and (iii) short wavelength negative energy 
modes, for which the geometric acoustics approximation 
may be valid, increase in number. The combination of 
these effects indicates that in the limit of high winding 
number it may be possible to recover a classical fluid de- 
scription of this system, for which the prediction of a 
thermal spectrum from Unruh [2] and Visser [18] should 
be experimentally verifiable. 

In contrast, for the relatively low winding numbers we 
have considered here, quantum effects are significant and 
the semiclassical approximation breaks down. It appar- 
ently becomes necessary to revise the concept of the ana- 
logue Hawking effect for trapped Bose-Einstein conden- 
sates in this regime. 
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